Energy landscape of a simple model for strong liquids 
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We calculate the statistical properties of the energy landscape of a minimal model for strong 
network-forming liquids. Dynamic and thermodynamic properties of this model can be computed 
with arbitrary precision even at low temperatures. A degenerate disordered ground state and log- 
arithmic statistics for the local minima energy distribution are the landscape signatures of strong 
liquid behavior. Differences from fragile liquid properties are attributed to the presence of a dis- 
crete energy scale, provided by the particle bonds, and to the intrinsic degeneracy of topologically 
disordered networks. 
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In recent years, the study of the statistical properties of 
the potential energy surface or landscape (PES) 0, Q 
sampled by liquids in supercooled states has received con- 
siderable attention, in the attempt to quantify the ther- 
modynamic properties of supercooled liquids and glasses 
QLIBIaul- The PES is analyzed to estimate the number 
and energy distribution Vl(E) of the local minima and 
the volume in configuration space sampled in vibrational 
motions around each local minimum. 

According to Angell's classification |8|, a liquid is 
named "fragile" if the temperature dependence of its 
transport coefficients shows large deviations from Arrhe- 
nius behavior. If no deviations are observed, the liquid 
is named "strong". For the case of fragile liquids — in 
the region of energies sampled during low temperature 
T equilibrium simulations with state of the art numer- 
ical resources — £l(E) is well described by a Gaussian 
distribution 0, 0, • The validity of Gaussian statistics 
in regions of the PES with energy lower than those ac- 
cessible in simulations, would suggest the existence of a 
finite Kauzmann temperature Tpc, where the configura- 
tional entro py S rnn f vanishes, and via the Adam-Gibbs 
relation [lfjllll| , a divergence of the characteristic relax- 
ation time at Tk- The quantification of the statistical 
properties of the PES for models of strong liquids is still 
under debate 0, It has been shown that, on low- 
ering T, deviations from Gaussian statistics take place. 
The configurational entropy does not seem to extrapolate 
to zero at a finite T, but the long equilibration times and 
the unknown value of the ground state energy prevent an 
unambiguous determination of the ground state degener- 
acy. A recent work of Heuer and coworkers fl3| suggests 
that the breakdown of Gaussian statistics originates from 
the emergence of a natural cut-off in fl(E), related to the 
formation of a fully connected network of bonds. 

In this Letter, we report a study aiming to clarify the 
statistical properties of the PES for a strong liquid, to 
provide a useful framework for interpreting results of real- 
istic models, and to shed light on the differences between 
fragile and strong liquids. We quantify the landscape 
properties for a simple model, similar in spirit to one 
previously introduced by Speedy and Debenedetti |l4j . 



We show that, with appropriate numerical techniques, 
both the dynamics and thermodynamic properties of this 
model can be computed with high accuracy and that no 
extrapolations are required to determine the low T be- 
havior. The strong liquid behavior is associated with the 
presence of two ingredients, which are characteristic of 
all network-forming liquids: a degenerate ground state 
and a discrete distribution of energies above the ground 
state, generated by the bond energy scale. 

We investigate a maximum valency model: a square 
well model of width A with a constraint on the maximum 
number of interacting particles. The interaction between 
two particles i, j that each have less than N max bonds to 
other particles, or between two particles already bonded 
to each other, is given by a square-well potential, 

{oo r < a 

-u Q a < r < A (1) 
r > <t + A. 

When i and/or j are already bonded to N max neighbors, 
then Vij(r) is simply a hard sphere (HS) interaction, 

10 r > a. 

The maximum number of bonds per particle is controlled 
by tuning N max . This model is particularly suited for 
theoretical and numerical studies. First of all, along an 
isochore, the system changes from a HS system (when 
ksT 3> Mo, ks is the Boltzmann constant) to a simple 
model for network- forming liquids (when fc^T < un) with 
coordination N max and no angular constraints |15I| . Sec- 
ond, the ground state energy for a system of N particles 
is known, being the energy of a fully connected network 
(—UoNN max /2). Third, as V(r) is simply a square well 
model, the local minima of the PES coincide with the 
bonding patterns of the system. Consequently, moving 
between minima takes place via breaking and reforming 
of bonds. Fourth, the volume in configuration space of 
each of these bonding patterns can be calculated with 
no approximation, as discussed in the following. In the 
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PES language, it means that the basin free energy Q can 
be calculated exactly. Finally, since the energy levels of 
the system are discrete by construction, the number of 
bonding patterns with M << NN max /2 broken bonds 
can be calculated by combinatorial factors. Exploiting 
these properties allows us to quantify precisely the sta- 
tistical properties of the PES of this model. 

We perform Monte Carlo (MC) and event driven 
molecular dynamics simulations by using A/(er + A) = 
0.03, Uo = 1, (7 = 1 and N max = 4. Entropy, S, is 
measured in units of ks- Setting ks = 1, energy E and 
temperature T are measured in units of uq. We study a 
system of N = 10 4 particles of equal mass m = 1, imple- 
menting periodic boundary conditions, at a fixed packing 
fraction = 0.30 0] from temperature T = 100 (where 
the HS limit is recovered) down to T = 0.04, where an 
almost fully connected network of bonds is obtained • 
No evidence of phase separation is observed [2(| • Fig. ^a, 
shows the T dependence of the potential energy per par- 
ticle E as & function of l/T. For a Gaussian PES [Hll22T|. 
E should be linear in l/T. Deviations from the l/T law 
are clearly seen at low T, when the energy approaches the 
ground state energy E gs = — 2, with a striking similarity 
with the T dependence of the energy of the sampled min- 
ima evaluated in landscape analyses of atomistic models 
for silica 0, ^| . The approach to E gs is well described 
by E — E gs ~ e~2T, providing an unambiguous way of 
calculating all thermodynamic properties down to T = 0. 

Next we provide evidence that the present model is 
a simple and satisfactory one for strong glass-forming 
liquids, by showing several features that are commonly 
observed in real systems. Fig. shows the T depen- 
dence of rj and D |23( . Both quantities display Arrhenius 
behavior at low T, as expected for strong liquids, with 
an activation energy controlled by the bond energy. The 
product -j? « , as expected from the Stokes-Einstein 
relation. We also evaluate the exponent (3 of the stretched 
exponential function that describes the long time decay of 
density autocorrelation functions, finding values (3 > 0.85 
for all the wavevectors. Such a large (3 value is in agree- 
ment with the experimentally observed direct correlation 
between (3 and strong behavior [24j . Finally, since at the 
investigated low T the potential energy has already ap- 
proached the ground state value (Fig. ^,), no significant 
drop in the specific heat will take place at the glass tran- 
sition temperature. This feature is often found in strong 
liquids [la.l25| . 

Next we evaluate the statistical properties of the PES. 
In the landscape approach, configuration space is parti- 
tioned into basins around the local minima of the PES. 
For the present model, different local minima can be un- 
ambiguously identified as different bonding patterns and 
the energy of the minimum coincides with the poten- 
tial energy of the configuration, i.e. with the number of 
bonds. The partition function can be written as a sum 
over all bonding patterns. At low T, it is convenient to 
write the sum over the number of M broken bonds, asso- 
ciating with each energy level E the number of distinct 
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FIG. 1: Temperature dependence of the potential energy per 
particle E (top) and of r\ (in units of (muo) 1 ^ 2 /a 2 ), D (in 
units of cr(wo/'7i) 1 ), and the product Dr//T (bottom). D is 
calculated from the long time limit of Sr 2 (t)/6t, with Sr 2 (t) 
the mean squared displacement, and rj as explained in |2ll |. 
In the top panel, the dashed line shows the l/T behavior 
expected for a Gaussian landscape [H l2^|. The inset shows 
the same data in a semi-log representation to highlight the 
validity of the low-T behavior hi[E — E gB ] = ln[2.56] — ^ 
indicated by the continuous line (see text). In the bottom 
panel the dashed lines are the best-fit Arrhenius behavior for 
r] and D. The continuous line corresponds to the expected 
value from the Stokes-Einstein relation (see text). 

network configurations with M = (E — E gs )N broken 
bonds (the degeneracy of the energy level) times the vol- 
ume in configuration space that can be sampled by the 
network without breaking or forming any bond (i.e. pre- 
serving the bonding pattern) . The logarithm of the num- 
ber of distinct networks with the same number of bonds 
defines S con f(E), while the logarithm of the sampled vol- 
ume defines the basin vibrational entropy S V ib ■ The sum 
Sconf + S V ib defines the total entropy of the model. 

In the present model, S V ib can be calculated without 
approximation by thermodynamic integration from a ref- 
erence Einstein crystal 26]. In this technique, a series 
of MC simulations are performed adding to the original 
Hamiltonian Ho a harmonic perturbation He{t N ]X) — 
A JT =1 (ri — f?) 2 , where (r?, ■■■,r%) is a typical configu- 
ration of energy E whose basin free energy needs to be 
evaluated, and the parameter A is the elastic constant of 
the harmonic perturbation. The system described by the 
Hamiltonian H(r N -,X) = H (f N ) + H E {r N ;\) is then 
simulated at fixed T and V with a MC technique, re- 
jecting all moves altering the bonding pattern. We study 
several values of A running from A = up to a value X max , 
whereupon the behavior of a system composed of 3./V in- 
dependent harmonic springs of elastic constant X max is 
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FIG. 2: Mean squared displacement per particle times A as 
a function of ln[\] for a typical basin explored at T = 0.04. 
The continuous line shows the harmonic behavior |fcsT. 

recovered. In this limit, the free energy is known ana- 
lytically. The vibrational free energy f v a (A = 0) of the 
basin to which r° belongs is given by 



fvib(T) = Fh e (T; Xmax) 

C^-^VnW, (3) 



where X^C^i - M ) i s the sum of the square displacements 
of all particles at a fixed value of A, brackets denote en- 
semble average, and Fh e (T; \ max ) is the free energy of 
the 3iV Einstein oscillators. Fig. [21 shows the A depen- 
dence of A(^ i (r i — f®) 2 )\/N for one specific state point. 
At large A values, the mean square displacement per par- 
ticle approaches the theoretically expected limit | . It 
must be noted that in the present model, since no energy 
changes are involved in sampling the configuration space 
volume at fixed bonding pattern, for each specific basin 
the basin free energy is purely entropic. 

Repeating the thermodynamic integration procedure 
for several different equilibrium configurations at differ- 
ent T, S V ib{E) can be precisely calculated. Fig. shows 
that the excess quantity, S%f b (E), over the ideal gas con- 
tribution, depends linearly on the basin depth E, a fea- 
ture shared with previously investigated models for su- 
percooled liquids 0, IE 



Syib(E) is well represented by 



the best-fit function: 



S% b = -8.8 + 4.27(E-E gs ). (4) 

Complementing these data with the calculation of the 
excess total entropy over the ideal gas contribution, S£% t , 
by integration of the T dependence of the specific heat 
over T from the known HS high T reference state [26) . 
allows us to evaluate S con f = S^ t — S^ b . Fig|3b shows 
that, close to the ground state, the E dependence of S con f 
is well described by a logarithmic combinatorial function: 



^conf 

(E) 



-2(E-E gs )ln[2{E-E : 



ys ) 



$conf \^gs ) 

- {E - E gs ). 



(5) 



The logarithmic term results from the number of dis- 
tinct ways (E — E gs )N bonds can be broken among 
N particles. A one-parameter fit (see Fig. HJJj) sug- 
gests S conf {E gs ) w 2.1 27]. Eqs. HE fully define the 




FIG. 3: Energy E [panels (a) and (b)] and T [panels (c) and 
(d)] dependence of Slot, S^ and Sconf - A typical error bar is 
shown in panel (b). Dashed lines for S%Z n f(E) and S^i b (E) are 
fits to Eqs l4l5l (see text). The continuous line for Slot is the 
result of summing the obtained fit functions. The continuous 
fines for Sfot(T), S^f b (T) and S con f(T) are calculated from 
the T dependence of E (see text). 



statistical properties of the PES. All thermodynamic 
functions can be evaluated from them. Indeed, solving 
d(S conf + Sl* b )/dE = 1/T, the T dependence of E is 
found to be E(T) ~ E gs = 2.56e~2r, providing the con- 
tinuous line in the inset of Fig. and parametrically, 
the continuous lines for S ex (T), S*f b {T) and S conf {T), 
compared with the numerical estimates in Figs. |3>|3Ji. 

According to the picture emerging from this study, 
strong liquid behavior is connected to the existence of 
an energy scale (provided by the bond energy) which is 
discrete and dominant as compared to the energetic con- 
tributions coming from non-bonded next-nearest neigh- 
bor interactions |29| . It is also intimately connected to 
the existence of a significantly degenerate ground state, 
favoring the formation of highly bonded states which can 
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still entropically rearrange to form different bonding pat- 
terns with the same energy [3(j • The presence of the bond 
energy scale j2j| also determines a distribution of energy 
levels above the disordered ground state which are log- 
arithmically spaced. These results help rationalize pre- 
vious landscape analysis of realistic models of network- 
forming liquids |12| and the recent observation by Heuer 
and coworkers that the breakdown of Gaussian landscape 
statistics is associated with the formation of a fully con- 
nected defect- free network p| . While in realistic models 
the ground state energy is not known, and the very long 
equilibration times needed prevent an unambiguous de- 
termination of the ground state degeneracy, both these 
quantities can be calculated for the present model. 

In this model, the discrete energy states provide a truly 
degenerate ground state. In real systems, small differ- 
ences in energy exist, e.g. resulting from second neighbor 
interactions, between configurations that in our model 
would correspond to the same energy state. Hence, in 
real systems the "ground state" energy is smeared out, 
singling out a unique lowest energy state. Consequently, 



for ksT values smaller than this spread in energy, S con f 
approaches zero. Therefore, our classical model is strictly 
applicable at temperatures where quantum effects are 
negligible and for values of &bT much larger than 
the spread in energy of the ground state. 

The simplicity of the model and the possibility of cal- 
culating accurately the relevant thermodynamic proper- 
ties at all T makes this model a relevant candidate for 
deepening our understanding of the differences between 
strong and fragile liquids. Results reported here suggest 
that strong and fragile liquids are characterized by signif- 
icant differences in their PES properties. In a simplified 
picture, a non-degenerate disordered ground state and 
Gaussian statistics characterize fragile liquids, while a de- 
generate disordered ground state and logarithmic statis- 
tics are associated with strong liquids. The origin of this 
difference is traced to the presence of a discrete energy 
scale, provided by particle bonds, and to the intrinsic 
degeneracy of topologically disordered networks. 
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